Validity of the scattering length approximation in strongly interacting Fermi systems 
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We investigate the energy spectrum of systems of two, three and four spin-^ fermions with short 
range attractive interactions both exactly, and within the scattering length approximation. The 
formation of molecular bound states and the ferromagnetic transition of the excited scattering state 
are examined systematically as a function of the 2-body scattering length. Identification of the 
upper branch (scattering states) is discussed and a general approach valid for systems with many 
particles is given. We show that an adiabatic ferromagnetic transition occurs, but at a critical 
transition point fci?a much higher than predicted from previous calculations, almost all of which 
use the scattering length approximation. In the 4-particle system the discrepancy is a factor of 2. 
The exact critical interaction strength calculated in the 4-particle system is consistent with that 
reported by experiment. To make comparisons with the adiabatic transition, we study the quench 
dynamics of the pairing instability using the eigenstate wavefunctions. 

PACS numbers: Valid PACS appear here 



I. INTRODUCTION 

The control afforded by Feshbach resonance phenom- 
ena in ultracold atomic gases has enabled the exploration 
of strongly correlated degenerate Fermi systems. In a re- 
cent study of the possibility of itinerant ferromagnetism 
Q-Q, Jo et al. [j| attempted to observe the physics be- 
hind the Stoner model in an atomic gas of 6 Li atoms. Ev- 
idence for ferromagnetic ordering was seen. This exper- 
iment has generated a great deal of theoretical research 
[3-fllj . The results have been debated as to whether a 
ferromagnetic transition or a strong correlation effect was 
seen. Quantitative comparison with experiment has not 
been achieved. 

A key issue is to hnd an appropriate model for the 
experiment. In the experiment a strong attractive inter- 
action is quickly switched on. Predictions of the criti- 
cal ratio of interaction strength to interatomic spacing 
for the ferromagnetic transition from mean field theory 
d, Hj], second order corrections [3, 0] and QMC calcu- 
lations [8 10] are on the order of kpa ~ 1; about two 
times lower than that from the Jo et al. experiment. In 
almost all calculations, a positive interaction Q or a 
Jastrow factor with two-body nodes [l(| is assumed, us- 
ing the scattering length approximation(SLA). Moreover, 
the various approaches differ in details of the nature of 
the transition p| H[ • 

The two-body SLA neglects the low-lying molecular 
states. The zero-energy s-wave scattering length a is de- 
fined by the long distance form of the out-going scattering 
wave 
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For a contact (zero range) potential, a is the radius of 
the first wavefunction node: ip (r = a) = 0. The SLA 



replaces the underlying atomic interaction by a purely 
repulsive potential which has the same two-body scatter- 
ing length. 

This is analogous to the idea of pseudo-potentials in 
electronic structure. A pseudo-potential can be gen- 
erated in an atomic calculation to replace the strong 
Coulomb potential of the nucleus and the effects of the 
tightly bound core electrons by an effective ionic poten- 
tial acting on the valence electrons and then used to com- 
pute properties of valence electrons in molecules or solids, 
since the core states remain almost unchanged. The ap- 
proach is widely used in electronic structure calculations. 
However, it leads to an inaccuate model if a pseudopo- 
tential is used for systems compressed to high density 
and the electron cores start to overlap. 

Many experiments in cold atomic systems are per- 
formed near Feshbach resonance where the scattering 
length is comparable to interatomic separation. In this 
situation, the lower-lying molecular bound states giv- 
ing rise to resonance can overlap, causing the scattering 
states to distort in order to remain orthogonal to the 
bound states. With more experimental effort expected 
in the study of related systems, precise and reliable com- 
parisons from quantum simulations will be important. 
Yet accurate many-body calculations will not be possi- 
ble without a quantitative understanding of the effective 
interactions and their effect on the different states. Even 
the identification of the scattering state in a dense system 
requires explanation. 

In this article, we quantify this effect by explicitly in- 
cluding the molecular bound states and treating the in- 
teraction exactly. We consider systems of two, three and 
four spin-i fermions interacting through a contact inter- 
action. The energy spectrum as a function of the two- 
body interaction strength is obtained by using an exact 
numerical method on a lattice and then extrapolated to 
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the continuum limit. We show how the upper branch can 
be identified for a many-body system. The properties of 
the nodal surface of the scattering many-body states are 
investigated. To compare with the exact solutions, cal- 
culations are also made with the SLA, by replacing the 
attractive contact interaction with a zero boundary con- 
dition. In both cases an adiabatic ferromagnetic tran- 
sition is stabilized. The SLA breaks down for large a, 
leading to a severe underestimation (by almost a factor 
of 2) of the transition point. 



II. METHOD 

We consider a system of two-component fermions mov- 
ing in a periodic box with length L to model a gas of 
6 Li atoms with two hyperfine species at non-zero density. 
All lengths are expressed in units of L and all energies 
in units of Kq = (^-) 2 - In the case a > ro (where 
ro is a measure of the effective interaction range), the 
interatomic potential can be modeled as a regularized 
(^-function: 
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where a is the zero-energy scattering length and m is the 
mass. We solve the Schrodinger equation by putting the 
system on a lattice with n points in each direction and 
recover the continuum limit by extrapolation. We then 
approximate the kinetic energy by two different discrete 
Laplacian operators [l3j: (1) the Hubbard model with 
nearest neighbor hopping and (2) a long range hopping 
model including up to the next nearest neighbors. We 
model the bare two-particle interaction by a point contact 
potential on the grid 
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where A = L/n is the grid spacing. Here U > is the 
strength of the attractive interaction; on the repulsive 
side of resonance U > Uoo, we have positive scattering 
length for unpaired atoms and the mapping relation be- 
tween the grid and continuum is [15 1 
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where the unitarity point a — > oo occurs at U^ 1 = 
(27r) -3 J" d 3 k(2e k ) _1 = jm/(h 2 A). Here e k is the single 
particle dispersion relation and 7 is a numerical constant 
defined by the discrete Laplacian. For choice (1) above, 
7 w 0.2527; for choice (2), 7 « 0.2190. When only near- 
est neighbor hopping is included, our Hamiltonian is the 
standard attractive Hubbard model, but scaled by 1/A 2 . 
Note our U value scales as A, while in the notation of 
the Hubbard model, is a constant. 

In the SLA, U has the opposite sign. In particular, 
when the scattering length a is large, Eq. ^ is replaced 




FIG. 1: (Color online) The two-body s-wave scattering energy 
Ei as a function of foa for grid sizes: n — 6,8, 10, 12, 20, 40. 
The inset shows the scaling with respect to l/n at kpa — 2.0 
with grid sizes up to n = 90. The long-range hopping model 
converges to the continuum limit faster than the Hubbard 
model. 



by a hard-sphere potential with radius a. If a is smaller 
than the grid spacing, a finite but negative value of U 
can be used in the SLA, leading to the repulsive Hubbard 
model, which clearly has a different strongly interacting 
or large a limit fll] from that of Eq. ([3]). 

To determine the eigenvalues and eigenstates, we start 
from a set of random trial states \ip^) where 1 < a < M, 
and evolve the states \i[> l+1 ) = (1 — tH)\iJj 1 ). At each 
step of the evolution, the state vectors are properly sym- 
metrized and orthogonalized. As i — > 00, the states con- 
verge to the lowest M eigenstates of the Hamiltonian 
H within a given symmetry sector. The errors are con- 
trolled and can be reduced arbitrarily with increasing the 
number of grid points or number of iterations. As dis- 
cussed below, the computational cost grows rapidly with 
system size, but significant reduction can be achieved by 
invoking symmetries. 

To assess the accuracy of iterative diagonalization, we 
test the method on a two-particle problem. The en- 
ergy of the lowest two-particle scattering state is plotted 
in Fig ([T]) as a function of the dimensionless scattering 
length kpa, where hp = {Sit 2 p) 1 ^ is the Fermi wave 
vector and p the particle density. Both discrete represen- 
tations of the kinetic energy operator were used and they 
converge to the same continuum limit: n —> 00; the long- 
range hopping is found to be less sensitive to the lattice 
spacing. Solving the two-particle problem also enables 
us to construct repulsive pseudo-potentials in the SLA 
by inverting the 2-particle Schrodinger equation. 
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FIG. 2: (Color online) Nodal surface for the scattering wave- 
function in a potential generated by two fixed particles located 
at (id/2,0,0) in infinite space with a/d = 1/10,1/3,1,5/2 
(expanding outward), where d denotes the distance between 
the two fixed scatterers. The solid (blue) lines correspond 
to the nodes in SLA and the dashed (red) lines to the exact 
nodes. SLA gives a reasonable approximation to the nodal 
surface for a/d < 1 but the deviation becomes significant for 
large scattering lengths. 

III. TWO FIXED POINT POTENTIALS 

The simplest case where the scattering length approxi- 
mation may fail is the scattering of a single particle off of 
two fixed contact potentials. This problem can be solved 
exactly in infinite space [llj]. The nodal surface of the 
zero-energy scattering state is given as the solution of: 



|r-Ri| |r-R 2 | a |Ri -R 2 |' v ' 

where Hi and R 2 are the location of the two fixed scat- 
terers. In the SLA, one would model the state by the 
ground state with nodes defined by |r — Ri| = a and 
|r — R 2 | = a. The nodal surfaces described by Eq. §5§ 
are shown in Fig ([2]) in comparison with corresponding 
spheres in SLA. Clearly the deviation from SLA becomes 
significant asa~ |Ri — R 2 |. In particular, the spherical 
surfaces in SLA becomes infinitely large at unitarity limit 
while Eq. §5§ gives rise to a finite surface. This result 
suggests that introducing a node in the two-body Jas- 
trow factor in the form f(r) ~ (1 — a/r) is insuffcient to 
characterize the effective pairwise repulsion on the upper- 
branch [lOj, as further discussed below in Sec I VII 

To study the effect of the SLA at finite density, we 
solved the same problem numerically in a finite periodic 
box. The results are summarized in Fig ([3]). It can be 
seen from the graph that at large scattering length (i.e. 
at high density) , the SLA significantly overestimates the 
scattering energy for the 3-body problem, i.e. the effect of 
low-lying molecular states cannot be ignored. The exact 
solution achieves a lower energy by distorting the nodal 



surfaces away from the union of two spheres required by 
the SLA. As we show below, this also applies to a system 
of more fcrmions. 



IV. FOUR- PARTICLE MODEL 

Now consider a minimal model for the ferromagnetic 
transition: four spin-i atoms in a cube with periodic 
boundary conditions and interacting with a contact po- 
tential. The spin polarized state #(1234) = ^(1234) ® 
Itttt) nas a totally antisymmetric spatial part ^(1234): 
for a contact interaction it is noninteracting with an en- 
ergy of 4.Ko in a zero total momentum eigenstate that 
has translational invariance. On the other hand, there 
are two linearly independent spin states with zero total 
spin S 2 = corresponding to unpolarized states: 

xms cx itm> + 144-tt) - \ [ in) + lit) ] ® [ in) + in) ] , 
xma oc[ it;) -in)]® [in) -in)]- 

The wavefunction is a linear combination of the above 
two states #(1234) = ip M A ® Xms + ipMS <8> xma- The 
symmetries of i(>ma and ipMS m coordinate-space are 
determined by the total antisymmetry of the complete 
wavefunction including spins and coordinates. 

For a system of four particles on a grid with n points 
in each direction, the discretized configuration space has 
n 12 grid points. Translational symmetries along the three 
spatial directions reduce the size of the configuration 
space by a factor of n 3 . Cubic symmetry of the peri- 
odic box reduces the number of independent wavefunc- 
tion values by a factor of 48 and permutation symmetries 
give an additional 24-fold reduction. We evolve pairs of 
non- magnetic states {4>mSj' 1 Pma} within the reduced do- 
main, and whenever the value of the wavefunctions on a 
grid point outside the reduced domain is needed in off- 
diagonal projections, the exterior point is mapped back 
into the reduced domain by symmetry transformations. 

The ferromagnetic transition is identified as the cross- 
ing between the lowest singlet scattering state and the 
fully ferromagnetic state. To investigate the effect of us- 
ing the SLA in multi-particle scattering process, the at- 
tractive contact interaction is replaced by a zero bound- 
ary condition and the resulting critical ferromagnetic 
density is estimated. 

Fig ((31) shows the energy spectrum of a four-particle 
system for n = 10 as a function of the coupling coefficient 
U. In this calculation, the lowest 30 states were followed. 
Note that we only considered states with the same sym- 
metries as the ground state, i.e. with even parity with 
respect to reflection in x, y or z. The resulting energy 
levels can be classified into three categories by their be- 
havior at strong coupling: two-molecule states, molecule- 
atom-atom states and four-atom scattering states. Level 
avoiding [17| can be observed between states belonging 
to different categories. 
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FIG. 3: (Color online) The scattering energy of a parti- 
cle moving in the potential generated by two fixed parti- 
cles located at x — L/4 and x — 3L/4, with y — L/2 and 
z = L/2. The SLA is obtained by replacing each poten- 
tial by a hard sphere (zero boundary condition) with the 
same scattering length. The SLA gives accurate energies in 
weakly-interacting limit (a/L < 0.2) but overestimates the 
scattering energy for a comparable with L. The inset shows 
the nodal region (\ipi\ < 10~ 4 ) for the scattering states with 
a/L = 0. 1 (black), 0.2(red), 0.4(blue). The surfaces become 
noticeably non-spherical for a/L > 0.1. 



strongly interacting regime because of the level avoiding 
phenomena. As shown in the inset of Fig ([J| , if the cou- 
pling coefficient U is tuned toward the resonance Uaa, 
it is energetically more favorable to jump through the 
successive avoided crossings. Thus, the change in the 
scattering energy due to an adiabatic tuning of the inter- 
action can then be determined by following the excited 
branch curve. It is drawn in bold in Fig Q. 

There is another way to identify the upper branch 
(scattering states) quantitatively by using the momen- 
tum distribution and the pairing order. First, consider 
the wavefunctions for the relative motion of two par- 
ticles interacting through a large scattering length of 
Eq. $Z§. The zero-energy scattering state in coordi- 
nate space ^(r) oc r _1 — a" 1 takes the form ^(k) oc 
4irk~ 2 — (27r) 3 a -1 (5(k) in momentum space and diverges 
at k = 0. By contrast, the bound state ip(r) oc r~ 1 e~' r / a 
takes the form ?/>(k) oc 4-7r(fc 2 + a~ 2 ) _1 in momentum 
space and remains finite at k = 0. The momentum distri- 
bution n(k) for scattering states is different from bound 
states at k = 0: scattering states have a larger fraction 
of particle occupation at k = 0. 

We also consider the pairing order defined by: 



.92 = n(^2S ritTj 



(6) 



for each state \i/j a }- The quantity gi measures double 
occupancy, and is related to the energy slope: 



V. IDENTIFICATION OF THE SCATTERING 
STATES 



The formation of molecular bound states is character- 
ized by the binding energy diverging linearly as U — > oc. 
In particular, the ground state wavefunction can be 
approximately written as ^o(13)^o(24) — Vo (14)^0 (23) 
where ip® is the two-body bound state, and the ground 
state energy is approximately twice the two-particle bind- 
ing energy. As seen in Fig. (3) the two-molecule states 
have an energy slope (dE/dU) about twice as large as 
the molecule-atom-atom states. As U — > oc molecules 
become tightly bound; their energy spacings can be un- 
derstood in terms of colliding molecules. For a lattice 
model, in contrast to a continuum model, the greater the 
internal binding energy, the greater the total mass of the 
molecule [181 ] . 

The scattering state of strongly repulsive atoms is an 
excited branch and all cold atom experiments performed 
in this regime are metastable. In Jo et al. experiment, 
the magnetic field ramp (~ 4.5ms) is much slower com- 
pared to the characteristic time scale of atomic collisions 
h/koTp ~ [is, and marches toward the resonance from 
the repulsive side a > 0. At low density or in the weakly 
interacting regime, the four-atom scattering state ap- 
proaches the noninteracting line 2Ko and the SLA is an 
accurate approximation. But there are some difficulties 
in defining the scattering state at high density or in the 
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For the scattering state, double occupancy decreases 
monotonically as the interaction strength is increased 
(see e.g. Ref. [ll|). Thus the scattering state in each lat- 
tice system is characterized by a vanishing energy slope 
as U — > oo, 



.92 -> 0, 



dE a 
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as can be seen in Fig (fj|. The pairing density is also 
related to the tail of the momentum distribution, which 
describes the short range physics. At large k, the momen- 
tum distribution takes the form n(k) —¥ C/k , where the 
coefficient C is called the contact in the Tan relations 
19]. In the continuum limit A«d, the contact C can 
be related to the energy slope, and hence gi, through the 
adiabatic sweep theorem -£^r = —-rzz^C, which in our 
system yields: 



.92 = L 



7 ' 



A 

4™ 



(9) 



where 7 is the numerical constant appearing in the defi- 
nition of Uao. For bound states this gives a finite gi and 
a finite energy slope. 

Thus, in addition to the momentum distribution at 
k = 0, we can identify the scattering state from the other 
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FIG. 4: (Color online) The energy spectrum of the lowest 30 s- 
states of four fermions with contact interactions on a 3D grid 
with n — 10. The energy levels can be classified into two- 
molecule states, molecule-atom-atom states and four-atom 
states. The ferromagnetic transition can be identified as the 
lowest scattering state (heavy dark line) crossing the hori- 
zontal ferromagnetic line around U/Uoo ~ 1.07. The inset 
shows an enlargement of the lowest scattering state and the 
associated avoided crossings. 



states by the magnitude of g^- scattering states have 
a smaller fraction of double occupation = r\,-. Note 
that the contact C measures the local density of pairs 
[2Cj . The momentum distribution n(k — 0) and the pair 
parameter 52 are plotted in Fig ([5]) as functions of the 
energy for kpa — 0.8 ~ 1.3. Scattering states are, by 
definition, in the range E/Kq > 2 and can be identified 
by the peaks of n(k = 0) and low values of <?2- 



VI. COMPARISON OF THE SLA AND EXACT 
RESULTS 

The ferromagnetic transition in the four-atom system 
occurs when the scattering state energy equals the non- 
interacting energy, AKq. For a n — 10 grid, the transi- 
tion occurs at U/Uqc ~ 1.07. Shown in Fig (|6]) is the 
energy of the four-particle unpolarized scattering state 
as a function of the scattering length kpa on grids with 
n = 6, 8, 10, 12 and their extrapolation to the continuum 
limit, n = 00. Avoided crossings between levels appear as 
kinks. The excited scattering state from the solution of 
the four-particle problem crosses the ferromagnetic line 
at kpa w 1.8, which is in remarkable agreement with the 
reported experimental value of kpa — 1.9 ± 0.2 

Also shown is the scattering energy using the SLA; this 
gives a ferromagnetic transition at kpa ~ 1.08 ~ 1.09 
for grid sizes n = 8, 10, 12, consistent with previous 
calculations [3, U, 04l3|- The earliest Fixed- node dif- 



FIG. 5: (Color online) The momentum distribution n(k = 
0) and the pairing parameter g2 of the lowest 30 s-states of 
four fermions with contact interactions on a 3D grid with 
n — 10 versus energy. Scattering states have, by definition, 
E/Ko > 2 and can be identified by the large magnitude of 
n(k = 0) and small magnitude of g2 compared to the other 
states. The peak structure at the scattering state diminishes 
as the interaction parameter kpa increases. 



fusion Monte Carlo calculations employed the repulsive 
Poschl- Teller potential (kpa « 0.86) [8J, hard spheres 
or repulsive soft spheres (kpa ~ 0.82) Q and included 
backflow effects (kpa ~ 0.96) [l(J. For attractive inter- 
actions modeled by spherical square wells or attractive 
Poschl- Teller potential, either variational Monte Carlo 
(kpa ~ 0.86) [|[ or fixed-node diffusion Monte Carlo 
10] (kpa ~ 0.89) calculates the upper-branch metastable 
state by imposing a nodal condition in the many-body 
wave function. The nodal condition ensures that the 
calculation consist of unbound fermionic atoms and no 
dimers or other bound molecules, by introducing a Jas- 
trow factor in the form of the scattering solution of the 
attractive potential corresponding to positive energy. As 
shown in Sec. IIII[ the nodal structure obtained this way 
deviates significantly from the true nodes in the upper 
branch. This explains why all these calculations gave re- 
sults similar to those from repulsive potentials, and all 
of them reproduced the predicted behavior of the mean 
field theory and second order corrections. 

The discrepancy between the critical values of kpa re- 
flects the limitations of perturbation theory in the regime 
of strong coupling. Compared to repulsive potentials, 
using Jastow factor with nodes and including backflow 
effects for attractive potentials improves the result by 
making nontrivial modifications to the nodal structure, 
but still gives answers not qualitatively different from the 
repulsive potential, and fails to reveal the inadequacy of 
the SLA. 

These observations suggest that lower-lying molecular 
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FIG. 6: (Color online) The four particle scattering energy as 
a function of the scattering length. The energy of the ferro- 
magnetic state is shown as the dashed horizontal line. The ex- 
trapolation to continuum is performed based on calculations 

st 

on grids with n = 6, 8, 10, 12, which exhibits the 1 -order 
linear scaling 1/iV to high accuracy. The transition to the 
ferromagnetic phase occurs at kpa « 1.8 while the SLA gives 
the transition at /cpa ~ 1.08. The inset shows the scaling 
with respect to 1/n near the transition point = 1.75 and 
1.8. 



states are responsible for delaying the formation of the 
ferromagnetic phase. However, calculations with more 
than 4 atoms are needed to determine finite size ef- 
fects. Such calculations are not feasible with the current 
method but might be possible with stochastic methods. 



VII. DYNAMICS OF FOUR-PARTICLE MODEL 

Because of the limited lifetime of the strongly in- 
teracting gas, however, the magnetic field ramp in ex- 
periment is not adiabatic and can lead to different 
explanations [U, [22| . A recent work [23| takes into ac- 
count the effect of atom loss by including a fictitious 
three-body term in the effective Hamiltonian of the Fermi 
gas and found that the critical interaction strength re- 
quired to stabilize the ferromagnetic state increases sig- 
nificantly. A full T-matrix analysis [2l| suggests that the 
pairing instability can prevail over the ferromagntic in- 
stability and the experimentally measured atom loss rate 
can be qualitatively explained in terms of the growth rate 
of the pairing order parameter after a quench. 

Thus, it is an interesting problem to study the dy- 
namics of the pairing instability after a quench using the 
wave functions obtained in this work. Since the contact 
C is identified as the integral over space of the expecta- 
tion value of a local operator that measures the density 
of pairs (20j . we characterize the pairing instability by 




t/X 

FIG. 7: (Color online) The first oscillation period in the evo- 
lution of g2(t) after a quench from a non- interacting state 
to an interaction strength Ufo,, where r = h/tF- The inset 
shows the nonmonotonic behavior of the maximum value of 
gi (t) in the first oscillation period as a function of the final 
interaction strength kpa. This calculation in done for four 
particles on a grid with n = 12. 

the count of double occupancy gi in Eq @. To study 
the dynamics of the pair formation, we choose the initial 
state to be the unpolarized four-particle ground-state in 
the noninteracting limit and expanded in the basis of 
the lowest 16 eigenstates with the final interaction after 
the quench. The time evolution is then evaluated using 
the eigenstate expansion \ip(t)) — J2 m c m,e~^ Emt \4'm) ■ 
The pairing density g%(t) — (V>(t)|<?2|^(*)) is used to 
characterize the pair formation and the atom loss into 
molecules. The evolution diagram of <?2 (t) is shown in 
Fig ([?])• The nonmonotonic dependence of the maximum 
value of g2(t) in the first oscillation period on the final 
scattering length kpa after the quench is in qualitative 
agreement with experiment and the full T-matrix theory 

Sa- 



vin. SUMMARY 

In summary, we have assessed the accuracy of the scat- 
tering length approximation at high density or strong in- 
teraction kpa > 1. It is demonstrated that if molecular 
states mix with excitations, non-magnetic states are sta- 
bilized. Identification of the upper branch in many-body 
calculations is discussed. The corresponding nodal struc- 
tures of the states are examined. The calculated critical 
interaction strength kpa for ferromagnetic transition is 
shown to be underestimated by a factor of two under the 
scattering length approximation. Although we solved the 
problem only for 4 particles, this minimal model suffices 
to show that ignoring the molecular states with the scat- 
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tering length approximation leads to inaccurate results in 
the strongly interacting regime. Hence it leads to severe 
errors in many-body calculations. That we get very good 
agreement with experimental estimates is encouraging 
but could be a result of cancellation of errors between the 
4 particle system and the thermodynamic limit. We in- 
vestigated the dynamics of pair formation. Non-monotic 
behavior of the pairing parameter J2i<j ^i^j * s observed 



as a function of the final interaction strength kpa after a 
quench. 
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